* ********************************************************************************************
* The Effects of Weather Shocks on Economic Activity: What are the Channels of Impact?
* Sebastian Acevedo, Mico Mrkaic, Natalija Novta, Evgenia Pugacheva, Petia Topalova
* With support from Gavin Asdorian, Olivia Ma, Jilun Xing and Yuan Zeng 
* Replication files for Table 1, column 6
* ********************************************************************************************
clear all
set more off
set matsize 11000

cd "C:\Users\..\replication" // set your working directory to where the replication folder is saved

* the Y variable that will be used in regressions: ngdprpc
global Y ngdprpc

* Horizons
global k 7

* ****************************************************************************************
* Our data - create extra variables
* ****************************************************************************************
use "data/country_dataset.dta", clear
drop if year > 2016

* year squared (for country-specific time trend squared)
gen sqyear = year ^ 2

* create variable for region-year (ry), encode it to be numeric, to be used as fixed effects
gen tempvar1 = wdi_region + " " + string(year) if wdi_region != ""
encode tempvar1, generate(ry)
drop tempvar1

* create country-specific time trend (_IifsXyea*) and time trend squared (_IifsXsqy*)
* dummies for region-year(_Iry*), country (_Iifscode*) and year (_Iyear*)
xi i.ifscode*year i.ifscode*sqyear i.ry i.year


* calculate GDP log difference at different horizons
gen ln_$Y = ln($Y)
by ifscode (year), sort: gen d0_ln_$Y = ln_$Y - l.ln_$Y
forval t = 1/$k {
	by ifscode (year), sort: gen d`t'_ln_$Y = f`t'.ln_$Y - l.ln_$Y
}


* generate squared terms of temperature and of precipitation
foreach X in UDEL_pw10 UDEL_pw90 UDEL_pw50 CRU_pw90 CRU_pw10 CRU_pw50 {
	gen tmp_`X'2 = tmp_`X' ^ 2
	gen precip_`X'2 = precip_`X' ^ 2
}

* ****************************************************************************************
* Variables where to save results
* ****************************************************************************************

gen period = -1 in 1

foreach X in t t2 p p2 t_ae t_em t_lic p_ae p_em p_lic{
	gen `X'_coef = 0 in 1
	gen `X'_se = 0 in 1
	gen `X'_pvalue = 0 in 1
}

gen nobs = 0 in 1
gen ncountry = 0 in 1
gen r2 = 0 in 1

* ****************************************************************************************
* Program
* ****************************************************************************************

cap program drop lpmregsq
program define lpmregsq
syntax,  yy(string) tt(string) pp(string) options(string)

forval t = 0/$k {
	* for horizons 0 and 1 we do not have forwards
	if `t' <= 1 {
		reg d`t'_ln_`yy' `tt' `tt'2 `pp' `pp'2 l.`tt' l.`tt'2 l.`pp' l.`pp'2 l.d0_ln_`yy' `options', cluster(ifscode)
	}
	
	* for horizons t = 2 to k we have forwards
	else if `t' > 1 {
		local s = `t' - 1
		reg d`t'_ln_`yy' `tt' `tt'2 `pp' `pp'2 l.`tt' l.`tt'2 l.`pp' l.`pp'2 f(1/`s').`tt' f(1/`s').`tt'2 f(1/`s').`pp' f(1/`s').`pp'2 l.d0_ln_`yy' `options', cluster(ifscode)
	}
	
	replace period = `t' in `=`t'+2'
	
	* temperature
	replace t_coef = _b[`tt'] in `=`t'+2'
	replace t_se = _se[`tt'] in `=`t'+2'
	replace t_pvalue = 2 * ttail(e(df_r),abs(_b[`tt']/_se[`tt'])) in `=`t'+2'	
	
	* temperature squared
	replace t2_coef = _b[`tt'2] in `=`t'+2'
	replace t2_se = _se[`tt'2] in `=`t'+2'
	replace t2_pvalue = 2 * ttail(e(df_r),abs(_b[`tt'2]/_se[`tt'2])) in `=`t'+2'	

	* precipitation
	replace p_coef = _b[`pp'] in `=`t'+2'
	replace p_se = _se[`pp'] in `=`t'+2'
	replace p_pvalue = 2 * ttail(e(df_r),abs(_b[`pp']/_se[`pp'])) in `=`t'+2'	
	
	* precipitation squared
	replace p2_coef = _b[`pp'2] in `=`t'+2'
	replace p2_se = _se[`pp'2] in `=`t'+2'
	replace p2_pvalue = 2 * ttail(e(df_r),abs(_b[`pp'2]/_se[`pp'2])) in `=`t'+2'	
	
	replace nobs = `e(N)' in `=`t'+2'
	replace ncountry = `e(N_clust)' in `=`t'+2'
	replace r2 = `e(r2_a)' in `=`t'+2'
	
	*** Temperature
	* lincom AE
	lincom _b[`tt'] + (2 * _b[`tt'2] * 11)
	replace t_ae_coef = r(estimate) in `=`t'+2'
	replace t_ae_se = r(se) in `=`t'+2'
	replace t_ae_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'

	* lincom EM
	lincom _b[`tt'] + (2 * _b[`tt'2] * 22)
	replace t_em_coef = r(estimate) in `=`t'+2'
	replace t_em_se = r(se) in `=`t'+2'
	replace t_em_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'

	* lincom LIC
	lincom _b[`tt'] + (2 * _b[`tt'2] * 25)
	replace t_lic_coef = r(estimate) in `=`t'+2'
	replace t_lic_se = r(se) in `=`t'+2'
	replace t_lic_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
	
	*** Precipitation
	* lincom AE
	lincom _b[`pp'] + (2 * _b[`pp'2] * 8)
	replace p_ae_coef = r(estimate) in `=`t'+2'
	replace p_ae_se = r(se) in `=`t'+2'
	replace p_ae_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'

	* lincom EM
	lincom _b[`pp'] + (2 * _b[`pp'2] * 9)
	replace p_em_coef = r(estimate) in `=`t'+2'
	replace p_em_se = r(se) in `=`t'+2'
	replace p_em_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'

	* lincom LIC
	lincom _b[`pp'] + (2 * _b[`pp'2] * 11)
	replace p_lic_coef = r(estimate) in `=`t'+2'
	replace p_lic_se = r(se) in `=`t'+2'
	replace p_lic_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
}

end


* ****************************************************************************************
* Run: rusults are stored in c
* ****************************************************************************************

* Column 6: Teulings and Zubanov (2014)
lpmregsq, yy($Y) tt(tmp_CRU_pw50) pp(precip_CRU_pw50) options("_Iifscode* _Iry*")
export excel period - r2 in 1/9 using "output/Table_1.xlsx", sheet("column_6") sheetreplace firstrow(var)
